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OBJECTIVES: Herniated discs and degenerative disc disease are major health problems worldwide. However, 
their pathogenesis remains obscure. This study aimed to explore the molecular mechanisms of these ailments 
and to identify underlying therapeutic targets. 

MATERIAL AND METHODS: Using the GSE23130 microarray datasets downloaded from the Gene Expression 
Omnibus database, differentially co-expressed genes and links were identified using the differentially co- 
expressed gene and link method with a false discovery rate <0.25 as a significant threshold. Subsequently, the 
underlying molecular mechanisms of the differential co-expression of these genes were investigated using 
Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis. In addition, the transcriptional 
regulatory relationship was also investigated. 

RESULTS: Through the analysis of the gene expression profiles of different specimens from patients with these 
diseases, 539 differentially co-expressed genes were identified for these ailments. The ten most significant 
signaling pathways involving the differentially co-expressed genes were identified by enrichment analysis. 
Among these pathways, apoptosis and extracellular matrix-receptor interaction pathways have been reported 
to be related to these diseases. A total of 62 pairs of regulatory relationships between transcription factors and 
their target genes were identified as critical for the pathogenesis of these diseases. 

CONCLUSION: The results of our study will help to identify the mechanisms responsible for herniated discs and 
degenerative disc disease and provides a theoretical basis for further therapeutic study. 
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■ INTRODUCTION 

The high prevalence of low back pain (LBP) makes disc 
degeneration a leading health problem. Currently, 70-80% of 
the adult population reports at least one episode of LBP 
during their life (1). In the American population, the 
incidence of symptomatic lumbar disc herniation is 1-2%, 
and approximately 200,000 lumbar discectomies are per- 
formed annually (2). Although there are several causes of 
spinal pain, intervertebral discs (IVDs) appear to be the 
most common source of chronic or recurring axial LBP (3). 
IVD degeneration is defined as an aberrant cell-mediated 
response to progressive structural failure and is a multi- 
factorial process influenced by genetics, lifestyle, and 
biomechanics (4). 
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Normal adult IVDs are made up of avascular tissue 
composed of a large amount of extracellular matrix (ECM) 
with a low cell density (approximately 1% of the total 
volume). Small solutes, such as glucose, lactic acid, and 
oxygen, in the inner parts of both the annulus fibrosus and 
nucleus pulposus are mainly transported by diffusion. The 
disc cells are responsible for maintaining the integrity of the 
ECM (5). Changes in the metabolic balance of IVD cells 
affect the quality and quantity of the ECM and its functional 
properties, and these changes can therefore be related to 
disc degeneration (6). Microscopically, disc degeneration is 
characterized by cell proliferation, cell necrosis or apoptosis, 
mucous degeneration, granular changes, and ingrowths of 
nerves or blood vessels (4,6). Cell proliferation leads to 
cluster formation, particularly in the nucleus pulposus (7). 
A recent study indicated that genetics plays a crucial role in 
disc degeneration, explaining 29-54% of the variance in 
adult populations (8). Environmental factors, such as 
mechanical loading and smoking, account for the remaining 
variance, although not all factors have been identified. 

In this study, the gene expression profiles of cells from 
herniated discs and discs with signs of degenerative disc 
disease, as well as data on transcription factors (TFs) and 
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target genes, were analyzed to identify underlying signaling 
pathways and critical molecules that regulate these path- 
ways at the transcriptional level. Recently, a number of 
methods, such as comparison and cluster analysis, have 
been adopted to analyze gene expression. However, many 
studies have reported that the identification of differentially 
expressed genes using these analyses cannot reveal the 
connections among various genes, and therefore, it is 
important to obtain clusters of genes that are differentially 
co-expressed, share common attributes and function colla- 
boratively (9,10). Using a variety of analytical methods, 
differentially co-expressed genes (DCGs) in cells from 
herniated discs and from discs exhibiting signs of degen- 
erative disc disease were analyzed using the DCGL package 
(11,12) in R. 

■ MATERIALS AND METHODS 

Microarray analysis 

The transcription profile of GSE23130 was obtained from 
the National Center for Biotechnology Information (NCBI) 
Gene Expression Omnibus (GEO, http:/ / www.ncbi.nlm.nih. 
gov/geo/). Based on the GPL1352 [U133_X3P] Affymetrix 
Human X3P Array, the RNA expression profiles of IVD cells 
from 23 patients with herniated discs or degenerative disc 
disease were analyzed. The 23 tissue samples for the 
microarray study consisted of one grade I, five grade II, nine 
grade III, five grade IV, and three grade V herniated discs or 
degenerative disc disease. All tissue samples were obtained 
from the Cancer Institute Cooperative Tissue Network 
(CHTN) or during surgeries. 

The Affymetrix GeneChip Operating System (GCOS) was 
used to calculate expression summaries with a target 
intensity of 100 using Microarray Suite version 5.0 (MAS 
5.0) (13). For completeness, the raw expression datasets 
were processed into expression estimates using the robust 
multiarray average method (14). The arrays were normal- 
ized using global scaling. 

Differentially co-expressed genes (DCGs) and links 
(DCLs) were identified using the DCGL method (http:// 
cran.r-project.org/ web /packages /DCGL/ index.html) (11,12), 
which was released as an R package that included two gene- 
filtering functions (an expression-based filter and a variance- 
based filter), three link-filtering functions (a systematic link 
filter, a percent link filter and a qlink filter) and five differential 
coexpression analysis functions (log ratio of connections 
[LRC], average specific connection [ASC], weighted gene 
coexpression network analysis [WGCNA], differential coex- 
pression profile [DCp], and differential coexpression enrich- 
ment [DCe]). These functions generally use gene expression 
matrices (with genes in rows and samples in columns) as the 
major input, and the ultimate output is genes ranked by p- 
value. The ^-values were adjusted for multiple comparisons 
using the false discovery rate (FDR) of Benjamini and 
Hochberg (15). An FDR-corrected p- values of 0.25 was used 
as the threshold value for screening the DCGs. The DCe has an 
additional output of classified DCLs. 

KEGG pathway enrichment analysis 

The Kyoto Encyclopedia of Genes and Genomes (KEGG) 
is a collection of online databases of genomes, enzymatic 
pathways, and biological chemicals (16). The PATHWAY 
database records networks of molecular interactions in the 
cells and variants that are specific to particular organisms 



(http://www.genome.jp/kegg/). A total of 130 pathways 
including 2,287 genes were collected from KEGG (updated 
on 2011.06). 

The Database for Annotation, Visualization and Integrated 
Discovery (DAVID) (17), a high- throughput and integrated 
data-mining environment, analyzes gene lists derived from 
high- throughput genomic experiments. DAVID was used to 
identify over-represented KEGG pathways based on the 
hypergeometric distribution, ^-values <0.05 were considered 
statistically significant. 

Analysis of transcription factors and their target 
genes 

According to the central dogma, a variety of methods can 
lead to differential gene expression. Gene expression is 
regulated by TFs at the transcriptional level. Therefore, the 
roles of TFs were further analyzed using the acquired data. 
Information on the chromosome region to which human TFs 
bind was obtained from the UCSC database (http:// 
genome.ucsc.edu). Subsequently, the chromosome annota- 
tions were downloaded from the NCBI database and 
analyzed to acquire the gene information of the above 
positions. Ultimately, a total of 210,000 regulatory pairs of 
TFs and their target genes were obtained. These pairs were 
then mapped to the above DCLs to determine the regulatory 
relationship between TFs and our DCGs. 

■ RESULTS 

Screening for differentially co-expressed genes 

Using the DCGL method described above, the gene 
expression profiles of cells from discs with signs of 
degenerative disc disease and from herniated discs with 
different degrees were analyzed, and genes with FDR<0.25 
were considered DCGs. A total of 539 DCGs and 113,866 
significantly related pairs of DCLs were selected. These 
findings suggest that these selected DCGs and their 
connections may have crucial roles in the development of 
herniated discs and degenerative disc disease. 

Analysis of biological pathways related to 
herniated discs and degenerative disc disease 

A total of 539 DCGs were included in the KEGG pathway 
enrichment analysis, and significantly altered pathways 
were ranked according to their p-values. The results showed 
that these DCGs were mainly involved in the following 
pathways (Table 1): protein processing in endoplasmic 
reticulum, Alzheimer's disease, focal adhesion, glycosphin- 
golipid biosynthesis-globo series, Huntington's disease, 
Parkinson's disease, oxidative phosphorylation, apoptosis, 
renal cell carcinoma, and selenoamino acid metabolism. 
Therefore, further analysis of these pathways may explain 
the pathogenesis of herniated discs and degenerative disc 
disease to some extent. The analysis of the apoptosis and 
ECM-receptor interaction pathways in particular may yield 
important insights because these pathways have been 
identified in previous studies. 

Analysis of gene regulatory relationships 

A total of 113,866 pairs of DCLs were matched with their 
already known TFs and their target genes, and 62 pairs of 
gene regulatory relationships (including 44 target genes) 
were found to be associated with herniated discs and 
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Table 1 - KEGG pathway enrichment analysis of 
differentially co-expressed genes. 



ID 



p-value Count Size 



Term 



4141 0.0011 11 167 Protein processing in endoplasmic reticulum 

5010 0.0012 11 168 Alzheimer's disease 

4510 0.0015 12 201 Focal adhesion 

603 0.0032 3 14 Glycosphingolipid biosynthesis - globo series 

5016 0.0074 10 184 Huntington's disease 

5012 0.0087 8 132 Parkinson's disease 

190 0.0094 8 134 Oxidative phosphorylation 

4210 0.0130 6 88 Apoptosis 

5211 0.0190 5 70 Renal cell carcinoma 

450 0.0191 3 26 Selenoamino acid metabolism 

5213 0.0275 4 52 Endometrial cancer 

4512 0.0380 5 84 E CM -receptor interaction 

5221 0.0433 4 60 Acute myeloid leukemia 

5215 0.0469 5 89 Prostate cancer 

511 0.0478 2 16 Other glycan degradation 

5210 0.0480 4 62 Colorectal cancer 

ID represents the identification number in the KEGG database; count 
represents the number of differentially co-expressed genes involved in 
this signaling pathway; size represents the number of total genes 
included in this pathway; and term represents the name of this KEGG 
pathway. 

degenerative disc disease (Figure 1). This result indicated 
that the presence of these 62 pairs of gene regulatory 
relationships may lead to the differential expression of 44 
target genes at the transcriptional level and ultimately result 
in the progression of disc herniation and of degenerative 
disc disease. 



■ DISCUSSION 

Based on our results described above, 539 DCGs were 
selected and shown to be enriched in several pathways in 
degenerative disc disease and herniated discs of different 
degrees. Among these pathways, apoptosis and ECM- 
receptor interaction (5) were demonstrated to be associated 
with these diseases. Both apoptosis and Fas/FasL expres- 
sion have been reported to be present in human disc 
specimens (18,19). The Fas/FasL system is a surface ligand 
that leads to rapid cell death when a Fas-expressing cell 
interacts with a second cell carrying the FasL receptor. Fas 
expression is not observed in normal discs but can be 
observed following disc injury, leading to disc degeneration 
(20). 

In addition, 113,866 pairs of DCLs were selected and 
matched with their known TFs and target genes. Finally, 62 
pairs of regulatory relationships that led to the identification 
of differentially expressed genes associated with degenera- 
tive disc disease and herniated discs of different degrees 
were obtained. There were numerous biological regulatory 
relationships between TFs and target genes (Figure 1). When 
the TFs were DCGs, some of their target genes were DCGs 
and some were not, indicating that many other factors were 
synergistically involved in the regulation of these target 
genes. However, although the TFs were not DCGs, all of 
their target genes were DCGs. These TFs may regulate their 
target genes through changes in their spatial structures and 
configurations and changes in localization. For example, 
NF-eB TFs are kept inactive in the cytoplasm by IeB 




Figure 1 - Networks of transcription factors and their target genes included in set of 539 differentially co-expressed genes. Blue dots 
represent target genes, and red dots represent transcription factors. The larger dots represent differentially co-expressed genes, and 
the smaller dots represent genes that are not differentially co-expressed. 
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(inhibitors of NF-eB) proteins (21). The activation process is 
regulated by the IeB kinase (IKK) complex. IKK activation 
leads to the phosphorylation and degradation of IeB 
proteins, subsequent NF-eB nuclear translocation and the 
transcriptional initiation of NF-eB-dependent genes (22). 
Although the amount of TF (NF-eB) is not changed, the 
activated TF can regulate its target genes and lead to their 
differential expression. 

In addition to these connections between TFs and target 
genes, many differentially expressed genes could still be 
regulated by TFs in indirect ways. For example, cathepsin K 
gene expression is significantly higher in more degenerated 
grade III to IV discs compared with healthier grade I to II 
discs (p = 0.001) (23). Cathepsin K, the protein encoded by 
CTSK, is a lysosomal cysteine protease involved in bone 
remodeling and resorption (24,25). Mice that are deficient in 
the cathepsin K gene have osteopetrosis of the long bones 
and vertebrae and abnormal joint morphology (26). 
Mutations in the human cathepsin K gene have been 
demonstrated to be associated with a rare skeletal dysplasia 
and pycnodysostosis (27,28). Cathepsin K over-expression 
could lead to Gaucher disease (29). In an analysis of CTSK, 
we found that CTSK could interact with the aryl hydro- 
carbon receptor (AHR) after interactions between eukaryotic 
translation initiation factor 2-alpha kinase 1 (EIF2AK) and 
ras homolog family member G (RHOG). The AHR protein 
contains several domains that are critical for its function and 
is classified as a member of the basic helix-loop -helix /Per- 
Arnt-Sim (bHLH/PAS) family of TFs (30,31). Therefore, 
although there is no direct relationship between CTSK and 
TFs, as shown in in Figure 1, CTSK can be regulated through 
interactions with other molecules (AHR- > RHOG — 
EIF2AK1— CTSK). It can be inferred that CTSK can be 
regulated by several types of upstream genes, resulting in 
differential expression. 

In conclusion, DCGs and the corresponding pathways, 
such as apoptosis and ECM-receptor interaction, were 
demonstrated to be involved in degenerative disc disease 
and herniated discs of different grades. Additionally, critical 
regulatory relationships between TFs and target genes, 
including one gene and several upstream regulators, were 
also identified in our study. Further experiments will be 
necessary to confirm these conclusions. 
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